#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <time.h>
#include <gmp.h> // Tem que colocar -lgmp na hora de compilar

#define N 25
#define BASE 10
#define MANTISSA 1024*1024*32

/*
Para um dual-core(2x1.8Ghz) com Mantissa de 4Kb:
	Um milhao de interacoes demora 101sec.

Mantissa de 8Kb:
	Um milhao de interacoes demora 273sec.

Mantissa de 16Kb:
	Um milhao de interacoes demora 753sec (12min).

Mantissa de 32Kb:
	Um milhao de interacoes demora 2127sec (~35min).

A complexidade da Mantissa é apx n*log(n)

Para 100mil interacoes com Mantissa de 32Kb:
	Demorou 210 segundos (~3.5 min)

Um aumento de 10x é linear. chute +9min - CONVERGE CERTO!
Mantissa de 64Kb 100mil interações: 549sec (~9.1min)

Mantissa de 64Kb:
	Um milhao de interacoes o chute eh 5700sec (~1:40h).


*/

void init_mpf_b(mpf_t rop){
	
	mpf_set_default_prec(MANTISSA);
	mpf_t x, y;

	mpf_init(x);
	mpf_init(y);

	// x=1
	mpf_set_ui(x, 1);
	// y=sqtr(2)
	mpf_sqrt_ui(y, 2);
	// x=1/sqrt(2)
	mpf_div(x, x, y);

	// return
	mpf_set(rop, x);

	mpf_clear(x);
	mpf_clear(y);
}

void media_simples(mpf_t rop, mpf_t a, mpf_t b){

	mpf_t aux;	
	mpf_init(aux);

	//aux=a+b;
	mpf_add(aux, a, b);
	//rop = aux/2
	mpf_div_ui(rop, aux, 2);

	mpf_clear(aux);
}

void media_geometrica(mpf_t rop, mpf_t a, mpf_t b){	

	mpf_t aux;
	mpf_init(aux);

	//aux=a+b;
	mpf_mul(aux, a, b);
	//rop = sqrt(aux)
	mpf_sqrt(rop, aux);

	mpf_clear(aux);
}

void equacao_t(mpf_t t, mpf_t x, mpf_t y, mpf_t a){
	
	mpf_t aux;
	mpf_init(aux);
	
	// t = t - x * pow((y-a),2);

	mpf_sub(aux, y, a);
	mpf_pow_ui(aux, aux, 2);
	mpf_mul(aux, aux, x);
	mpf_sub(t, t, aux);

	mpf_clear(aux);
}

void equacao_pi(mpf_t pi, mpf_t a, mpf_t b, mpf_t t){

	mpf_t aux;
	mpf_init(aux);
	
	//pi = pow((a+b),2)/(4.0*t);

	mpf_add(aux, a, b);
	mpf_pow_ui(aux, aux, 2);
	mpf_div(aux, aux, t);
	mpf_div_ui(pi, aux, 4);

	mpf_clear(aux);
}
int main(int argc, char* argv[]){
	
	unsigned long int i;
	time_t t_zero, t_final;
	double delta_t;
	
	// Marca o tempo
	time(&t_zero);
	
	// Float com precisao arbitraria
	mpf_set_default_prec(MANTISSA);
	mpf_t a, b, t, x, y, pi;

	mpf_init(a);
	mpf_init(b);
	mpf_init(t);
	mpf_init(x);
	mpf_init(y);
	mpf_init(pi);

	//	a = 1.0;
	//	b = (1.0/(sqrt(2)));
	//	t = (1.0/4.0);
	//	x = 1.0;
	//	pi=3.14, y=1.0; // Independnetes

	mpf_set_d(t, (1.0/4.0));
	mpf_set_ui(a, 1);
	mpf_set_ui(x, 1);
	init_mpf_b(b);

	/*
	 * 	A(n+1) = A(n)+B(n)/2
	 * 
	 *  B(n+1) = sqrt(A(n)+B(n))
	 * 
	 *	T(n+1) = T(n) - P(n)(A(n) - A(n+1))^2
	 *  
	 * 	P(n+1) = 2*P(n)
	 * 
	 * 	pi = ((A(n)+B(n))^2)/4T(n)
	 * 
	 * */

	for(i=0; i<N; i++){
		
			mpf_set(y, a);				//y = a;
			media_simples(a, a, b);  	//a = (a+b)/2.0;
			media_geometrica(b, b, y);	//b = sqrt(b*y);
			equacao_t(t, x, y, a);		//t = t - x * pow((y-a),2);
			mpf_mul_ui(x, x, 2);		//x *= 2;
			equacao_pi(pi, a, b, t);	//pi = pow((a+b),2)/(4.0*t);
	}

	// Limpa a Caquinha
	mpf_clear(a);
	mpf_clear(b);
	mpf_clear(y);
	mpf_clear(t);

	mpf_out_str(NULL, BASE, 0, pi);

	mpf_clear(pi);
	
	time(&t_final);
	delta_t = difftime(t_final, t_zero);
	printf("\nTempo = %.2lf segundos\n", delta_t);

	return 0;
}
